Synchrotron‐source micro‐x‐ray computed tomography for examining butterfly eyes

Abstract Comparative anatomy is an important tool for investigating evolutionary relationships among species, but the lack of scalable imaging tools and stains for rapidly mapping the microscale anatomies of related species poses a major impediment to using comparative anatomy approaches for identifying evolutionary adaptations. We describe a method using synchrotron source micro‐x‐ray computed tomography (syn‐μXCT) combined with machine learning algorithms for high‐throughput imaging of Lepidoptera (i.e., butterfly and moth) eyes. Our pipeline allows for imaging at rates of ~15 min/mm3 at 600 nm3 resolution. Image contrast is generated using standard electron microscopy labeling approaches (e.g., osmium tetroxide) that unbiasedly labels all cellular membranes in a species‐independent manner thus removing any barrier to imaging any species of interest. To demonstrate the power of the method, we analyzed the 3D morphologies of butterfly crystalline cones, a part of the visual system associated with acuity and sensitivity and found significant variation within six butterfly individuals. Despite this variation, a classic measure of optimization, the ratio of interommatidial angle to resolving power of ommatidia, largely agrees with early work on eye geometry across species. We show that this method can successfully be used to determine compound eye organization and crystalline cone morphology. Our novel pipeline provides for fast, scalable visualization and analysis of eye anatomies that can be applied to any arthropod species, enabling new questions about evolutionary adaptations of compound eyes and beyond.


| INTRODUC TI ON
There is a rich history of using insects to understand behavioral and anatomical diversity (Chown & Terblanche, 2006;Price et al., 2011).
Insects represent the largest group in the animal kingdom and their absolute numbers are also matched by their diversity in phenotypes, behavior, and anatomy (Stork, 2018).Classically, morphological variation that could be observed by the naked eye provided the necessary evidence for fundamental theories in evolution including natural selection, speciation, mimicry, and mate preference (Butler, 1963;Darwin, 1859;Poulton, 1909), to name a few.
More recently, the revolution in genetics and genomics has allowed for identifying genetic variation that drives variation in these observable traits (Baxter et al., 2010;Dobzhansky, 1982;Kronforst et al., 2006).However, microscopic studies have lagged behind, largely due to a lack of experimental tools to rapidly visualize and analyze fine structural detail over large volumes and algorithmic tools to analyze the resulting large image data sets with minimal human effort.While there has been a recent push to test different techniques for studying morphology, most methods do not provide a satisfactory balance between higher resolution and lower computational power (Friedrich et al., 2014;Van de Kamp et al., 2014;Wipfler et al., 2016).
Electron microscopy (EM) can provide the requisite resolution but is typically limited to scanned EM, (SEM) which visualizes external morphologies (Hao et al., 2023;Schwarz et al., 2011).A full 3D EM reconstruction using serial block face SEM, focused ion beam SEM, or transmission EM remains time-and computation-intensive.
We, and others, have recently shown that the sample preparation for EM using osmium tetroxide, which is species independent, provides excellent contrast in X-ray tomography microscopes (Dyer et al., 2017;Johnson et al., 2006;Ribi et al., 2008;Van den Boogert et al., 2018).Using X-ray tomography, large volumes of brains (even entire mouse brains) can be imaged in 3D at submicron resolution quickly (imaging rates of 0.067 mm 3 /min; Foxley et al., 2021).Here we demonstrate a pipeline for synchrotron source X-ray computed tomography (syn-μCT) performed at the Advanced Photon Source (APS) at Argonne National Laboratory (ANL) for high-throughput 3D imaging of the brains and intact eyes of a variety of butterflies.
2. We developed a novel embedding method that allows for automatically imaging multiple species eyes in a single imaging run to enable high-throughput imaging.
3. We developed a machine vision pipeline to extract the relevant morphological features from X-ray datasets and used these reconstructions to better understand microscopic variability in the morphology of cells in the light path across species.
4. Specifically, we analyze these new data sets in the context of pioneering work in Hymenoptera species (e.g., bees and parasitic wasps) that determined an optimal ratio of interommatidial angle to resolving power (Barlow, 1952).This ratio of interommatidial angle to resolving power is hereafter referred to as the "Barlow ratio" and is dimensionless as both angle and resolving power are in degrees.We extend this work by showing the Barlow ratio of the ommatidia in butterfly species falls near the theoretical optimum.By leveraging the full 3D datasets, we, however, find significant variation across an individual eye.
5. Finally, we use an amalgamation of individual crystalline cone measurements across individual eyes to generate a representative 3D crystalline cone for each sample within and across species.Generating the morphology of these cones allows for the mapping of light as it travels through this structure to the rhabdom.We observe cone shapes that vary both across the eye of an individual and between individuals (Figure 5).This technique allows for the dissection of these effects at fine detail across the eye and could support studies of cone optics and variation in and between species.

| MATERIAL S AND ME THODS
Samples from seven animals across six species of butterflies (Heliconius cydno, Strymon melinus, Calycopis cecrops, Polygonia interrogationis, Polites peckius, and two Pieris rapae) were prepared for electron microscopy (Hua et al., 2015) and assembled in plastic pillars vertically to stabilize the samples for imaging (Figure 1a (APS) at Argonne National Laboratory using syn-μCT using an automated z-axis tiling approach for unassisted imaging of multiple insect eyes (Figure 1b).The resulting X-ray data sets, with a total volume of 14.3 mm 3 and an isotropic resolution of ~0.6 μm resolved fine structure in the eye across all species, most notably the crystalline cones (Figure 1c, Figure S1).We next developed our analysis pipeline by focusing on the crystalline cones due to its notable variability across species upon visual inspection.
We used an analysis pipeline to extract the relevant features from the X-ray datasets.For example, Figure 2a shows the segmentation output of ilastik (Berg et al., 2019), a free open-source software for image classification and segmentation.The output from ilastik gave us clusters of points corresponding to each crystalline cone, which we analyzed in Matlab (MATLAB, 2021) and Python (Van Rossum & Drake, 2009; Figure 2b, Figure S1).

| Sample preparation
Insect samples were either collected in the wild in Chicago, IL (Pieris rapae, Polites peckius, and Polygonia interrogationis), collected from our breeding colonies at The University of Chicago (Heliconius cydno), or provided by Erica Westerman (University of Arkansas) (Strymon melinus and Calycopis cecrops).For dissections, insects were anesthetized by placing them at 4°C for ~10 min.Insects were then submerged in ice-cold phosphate-buffered saline (PBS) and dissected in PBS under a stereomicroscope to remove the cuticle outer layer and expose the brain.Brains with eyes intact were then cut from the body and submerged in fixative solution consisting of 0.1 M Sodium Cacodylate buffer, pH 7.4, 2% paraformaldehyde, and 2.5% glutaraldehyde.Brains were incubated in fixative for ~24 h, gently rocking at 4°C.The next day, brains with eyes were prepared using electron microscopy protocols as previously described (Hua et al., 2015).Briefly, brains were washed extensively in cacodylate buffer at room temperature and stained sequentially with 2% osmium tetroxide (EMS) in cacodylate buffer, 2.5% potassium ferrocyanide (Sigma-Aldrich), thiocarbohydrazide, unbuffered 2% osmium tetroxide, 1% uranyl acetate, and 0.66% Aspartic acid buffered Lead (II) Nitrate with extensive rinses between each step with the exception of potassium ferrocyanide.The samples were then dehydrated in ethanol and propylene oxide and infiltrated with 812 Epon resin (EMS, Mixture: 49% Embed 812, 28% DDSA, 21% NMA, and 2.0% DMP 30).Samples were cured in custom cylindrical molds to stack multiple brains into one sample and to remove any edges to the resin that may affect X-ray imaging.The resin-infiltrated tissue was cured at 60°C for 3 days.

| μ X-ray computed tomography
The syn-μCT data were acquired as previously described (Foxley et al., 2021).Briefly, we used the 32-ID beamline at the Advanced Photon Source, Argonne National Laboratory.The setup consists of a 1.8 cm-period undulator operated at a low deflection parameter value of K = 0.26.This yields a single quasi-monochromatic peak of energy 25 keV without the losses incurred by use of a crystal monochromator.For a sample 68 m from the undulator, this produces a photon fluence rate of about 1.8 × 107 photons s −1 μm −2 .The x-rays were imaged using a 10 μm thick thin-film LuAG:Ce scintillator producing visible-light images then magnified using a 10× Mitutoyo long working distance microscope objective onto a 1920 × 1200 pixel CMOS camera (Point Gray GS3-U3-51S5M-C).The effective object space pixel size was 600 nm isotropic.The thickness of the thin-film scintillator matched the depth of focus of the objective lens, achieving a spatial resolution equivalent to the resolving power of the lens (1.3 μm for a NA of 0.21).Since the camera field of view was substantially smaller than the sample, a mosaic strategy was employed (Vescovi et al., 2018).
The sample was mounted on an air-bearing rotary stage (PI-Micos UPR-160 AIR) with motorized x/y translation stages located underneath and x/y piezo stages on top.Typical exposure time for a single projection image at one mosaic grid point and one rotation angle was 30 ms. 360° rotation angles were used at each grid point.
The sample was translated through a 6 × 18 tomosaic grid.

| Data analysis
Crystalline cones from the raw x-ray datasets were segmented using the software ilastik and code based off cc3d (Silversmith, 2021).This generated sets of voxels corresponding to each of the crystalline cones.Outliers in the set of points that were not part of the cones were deleted manually.
We defined the center of each crystalline cone as its center of mass.Then we estimated the local radius of the eye by fitting a sphere to clusters of 60 points corresponding to the crystalline cone centers.We chose 60 points because this encompasses a hexagonal array surrounding a single point extending 4 ommatidia out in all directions.Vectors from the center of the sphere to the center of each cone were calculated.Once we have defined the 'center' of the eye from the local curvature we can then use the vectors from that putative center to the centers of the cones to define an ommatidial angle.The angles between a cone's vector and its six nearest neighbors' vectors were averaged, and this was used as the (local) interommatidial angle (Δɸ).The average distance to the six nearest neighbors was used as the diameter of the ommatidium (D).Since the center of each cone lies below the surface of each eye facet, this systematically underestimates the value of D by potentially a significant fraction of the cone length times Δɸ (measured in radians).This systematic error is then of order 2 μm or less, which is considerably smaller than both the mean and the variance of D (Table S1).
Resolving power was calculated by θ = 1.22*λ/Dwhere λ is the wavelength of light.For our analysis, λ = 500 nm, as it corresponds to broad peaks in both the typical sunlight spectrum and photoreceptor sensitivity in many insects.This is also the λ that Barlow used for his calculations.The ratio Δɸ/θ, aka the Barlow ratio, was also calculated.Extreme outliers were cut off when we noted corresponding defects in the x-ray images or where the values seemed biologically implausible (e.g., interommatidial angles greater than 90 degrees).These outliers occurred almost exclusively at the edges of the eyes.Cone shapes were determined by centering and overlaying all cones within a single eye, and keeping the collection of shared points, points with at least 25% overlap, across all cones.This was done to reduce noise in the segmentation of individual cones.The number of cones overlayed per eye varied from about 600 to 3000.
There appeared to be some variation in cones across the eye, but the biggest deviations from the average seemed to come from the cones at the outer edges of the eye.Then, the boundary of the cone was calculated from this set of overlapping points using the "boundary" function on Matlab.We used the Pearson correlation coefficient to determine the relationship between wingspan and cone length as well as wingspan and cone ratio.We also calculated the aspect ratio and cone ratio for each individual cone.

| RE SULTS
In order to understand whether ommatidial diameter, interommatidial angles, resolving power, and Barlow ratios changed within individuals, we first looked at how these different parameters varied across the surface of a single eye.When creating these maps, we found that ommatidial diameter (as well as resolving power) appeared to vary gradually from areas with larger diameters up to 35.3 μm to those with smaller diameters down to 17.1 μm (Figure 2c,d).In contrast, the changes in interommatidial angle and Barlow ratio across the eye were not as smooth, with transitions from higher to lower acuity areas being less clear by visual inspection (Figure 2e,f).
Next, we asked how these parameters varied across different individuals by reporting the statistics of the distributions for these variables for each individual.The average median ommatidial diameter across all individuals was 24.5 μm, ranging from 20.54 to 31.09 μm, with an average interquartile range of 3.36 (Figure 3a).
Resolving powers had an average median value of 1.46 degrees, ranging from 1.12 to 1.70 degrees, with an average interquartile range of 0.18 (Figure 3b).The median interommatidial angle measured across all species ranged from 1.42 to 1.87 degrees, with an average of 1.66 degrees and an average interquartile range of 1.01 degrees (Figure 3b).Medians for the Barlow ratio ranged from 0.846 to 1.49 with two individuals having medians within the optimal range.The average median Barlow ratio was 1.17 with an average interquartile range of 0.80 (Figure 4a).
Because our method provided great enough resolution to clearly distinguish whole crystalline cones, we decided to look at the micron scale morphology of this structure, which guides the light focused by the cornea and lens to the ommatidia.When looking at the size and shape of this structure, we found variation across individuals that would need to be disentangled from variation across the eye with larger datasets.For instance, the average "typical" cone length across species (Figure 5) was 44.4 μm, but ranged from 22.2 to 72.0 μm, and some cones had a defined point at the bottom while others were more rounded.In order to quantify how tapered a cone was, we took the ratio of the cone diameter at 10% of the length from the bottom and the maximum diameter at the top and found the mean of this ratio was 0.3669.
The most tapered cone had a ratio of 0.2121, while the least tapered cone was less than half as tapered with a ratio of 0.5362.
When analyzing morphological data, it is important to consider scaling relationships in the data (Jablonski et al., 1996).In our data, when examining allometric relationships, we found that there was a negative correlation between the wingspan of a species and the typical cone length (r(7) = −.8881,p = .0076)as well as a negative correlation between wingspan and cone ratio (r(7) = −.7949,p = .0326).Besides looking at the typical crystalline cones, we also looked at the shapes of all the cones across each eye.Changes in different parameters such as the length (Figure 5b), aspect ratio (Figure 5c), and cone ratio (Figure 5d) could be seen across the eyes of all individuals.

| DISCUSS ION
Our method provides a new way to study insect morphology, especially the individual components of the eye, using a higher contrast staining method, higher resolution syn-μCT, and a novel analysis pipeline.With the method, an entire eye can be surveyed for microscopic features like ommatidial diameters, angles, and cone morphologies.Previous methods measuring microscopic features have either imaged smaller volumes at higher resolution (e.g.Hao et al., 2023) or larger volumes at lower resolution (e.g.Currea et al., 2023).Using our approach, we were able to analyze seven insects of six different butterfly species to show that the Barlow ratio of the ommatidia falls in or near the theoretical optimum, but notable portions of the eye have Barlow ratios greater than this optimum.This suggests that portions of the visual scene are undersampled.In previous work that has shown a similar kind of spatial undersampling in insect eyes, it has been suggested that this is due to motion blur from the animal moving about its environment (Land, 1997).This indicates that theoretical models must also account for the angular velocity of the organism or objects in its visual environment and other processing that happens later in the visual system to optimize an insect's vision.Additional research is needed to assess how different luminance may change the Barlow ratio and if the theoretical models that account for light levels, as described by Snyder et al. (1977), are correct.

| Limitations
There are several limitations to this study.One issue is that our imaging often did not cover the entire eye.While this is not ideal, we were still able to analyze large enough portions of the eyes to capture the variability across eyes.This study is complementary to previous work that allows sampling of different parts of the eye.However, there is a future planned upgrade to the APS synchrotron that will enable imaging of entire eyes and nervous systems of insects (Argonne National Laboratory, n.d.).
Another limitation in study is that sample preparation for electron microscopy is well known to change the native structures of brain tissue (Zhang et al., 2017).However, most of these artifacts involve changes in the volume of the extracellular space (Pallotto et al., 2015;Van Harreveld & Steiner, 1970).We analyzed crystalline cones, which are composed of concentrated, hydrophobic proteins in closely related moths and likely less susceptible to dehydrationbased distortions (Schlamp, 1989).We designed an analysis pipeline robust to small changes in orientation, thereby preserving local curvature and diameters.Finally, we see smooth variation across individual eyes, which gives us confidence that the differences observed are not simply noise from artifacts.

| Comparison to prior work
Previous analyses of ommatidial diameter and interommatidial angles were done using light microscopy (and more recently fluorescence microscopy) and analyzed manually (Baumgärtner, 1928;del Portillo, 1936;Horridge, 1978;Rigosi et al., 2021;Rutowski & Warrant, 2002), and therefore would take a much longer time to collect data on the same volume of eye.Previous methods for calculating interommatidial angles include observing how many ommatidia pseudopupils crossed while rotating the eye a certain angle, using the optomotor response, and manually measuring histological sections (Baumgärtner, 1928;del Portillo, 1936;Götz, 1965;Horridge, 1978;Rigosi et al., 2021;Rutowski & Warrant, 2002).All of these methods are subject to human error, but our method provides an automated way to calculate both the interommatidial angle and the ommatidial diameter.Several newer methods have been proposed for measuring interommatidial angles and other eye parameters.One such method involves staining photoreceptors with fluorescent dyes to measure interommatidial angles using the pseudopupil in insects with dark eyes (Rigosi et al., 2021).One advantage of this fluorescence method is that it can be done using live animals and avoids any distortion that may occur during sample preparation.However, this is the only parameter that can be measured with this technique and cannot reveal the morphology of internal structures.
The μCT method has recently been used to measure angles and other eye parameters in bees (Taylor et al., 2019) and ommatidial diameters in other compound eyes (Currea et al., 2023), but our method using the 32-ID beamline achieves ~18× or ~ 170× greater resolution respectively, and imaging speeds of ~1 mm 3 /30 min.This enhanced resolution combined with our novel embedding method allows for greater automated throughput.Furthermore, conventional lab-based μCT imaging that can achieve comparable spatial resolution (Alba-Tercedor et al., 2021) have worse contrast resolution than syn-μCT (Goyens et al., 2018).Additionally, our staining method provides even greater contrast and allows us to better see crystalline cones, whereas previous μCT reconstructions were unable to capture this structure.We found large variation in the sizes and shapes of the typical crystalline cones across individuals and especially species.Using our 3D models of these cones, further research can be done to explore how light passes through these structures and impinges on the rhabdom.

| The Barlow ratio, variation, and motion blur
There was considerable variation in all measurements across individual eyes.All the species had average Barlow ratios near the theoretical optimum, but large portions of each eye had ratios that were greater than expected, meaning the interommatidial angle was greater than the resolving power, suggesting the visual scene is undersampled.Undersampling of the visual scene has been observed in other insects.In previous studies that have measured a similar ratio, the acceptance angle to the interommatidial angle, in other diurnal insects also found that the visual scene was undersampled (Land, 1997).One reason insect vision might be undersampled is to account for motion blur.For example, according to a paper from Snyder et al, the fly Musca would have a Barlow ratio of 2.13, which is greater than both Barlow's optimum and the optimum calculated by Snyder et al. (1977).However, this value did approach a value that Snyder et al. deemed more reasonable once angular velocity of the insect was accounted for.Finally, Snyder et al. also looked at how different light levels would affect the theoretical optimum for p = D*Δɸ = 0.61*(Δɸ/θ).They theorized that in lower light conditions, the optimal p would be larger.Further research must be done to examine crepuscular and nocturnal Lepidoptera to determine if this is indeed the case.
Finally, we determined the morphology of a typical crystalline cone for each species.We found considerable variation in height and width within and across individuals, which could be due to scaling with total body size or cone density within the eye.
Previous analyses have identified that the point-like end of the crystalline cone corresponds with the focal point of the lens, allowing the most efficient transfer of photons into a single rhabdom (Schwarz et al., 2011).For species adapted to low light, the hypothesis is instead that cones are larger and more bulbous with the focal point well inside the cone, which is believed to confer an advantage for greater light collection by transmission through a 'clear zone' to multiple rhabdoms (Warrant, 2017).Since the crystalline cone's function is to funnel light onto the rhabdom, further studies could potentially determine how cone optics vary across the eye and between species.Measurements of body size that incorporate forewing length are correlated with larger eye sizes (Seymoure et al., 2015), and longer cones correspond to smaller wingspans, suggesting smaller Lepidoptera have flatter lenses as they have a longer focal length.Syn-μCT also clearly shows the shape of the lens, so our method would be useful in testing this hypothesis.Future work will explore these differences more fully, by modeling the wave optics of light passing through lenses and cones with these different shapes.

CO N FLI C T O F I NTER E S T S TATEM ENT
We have no conflicts of interest to disclose.
), and large sections of eyes were imaged at the Advanced Photon Source F I G U R E 1 X-ray analysis pipeline showing (a) diagram of insect eyes stacked in a vertical column, (b) a diagram of the sample rotating and moving vertically in the X-ray beam, and (c) a raw X-ray image.Butterfly in (a) is from (Gallice, 2012).

F
I G U R E 2 (a) Segmented out crystalline cones with the inlay showing that cones are labeled as separate objects.(b) Centers of crystalline cones plotted in Matlab where the rest of our analysis took place.(c-f) show scatter plots showing how (c) ommatidial diameter, (d) resolving power, (e) interommatidial angle, and (f) Barlow ratios change across the eye.The portion of the eye shown here is from Polites peckius.Each point represents one ommatidium.

F
I G U R E 3 Plots showing the median (red line) (a) ommatidial diameter in micrometers, (b) resolving powers, (c) and interommatidial angles in degrees for the seven individuals.Boxes show interquartile range and whiskers show the lower and upper quartiles (Outliers are shown in Figure S2 for clarity of presentation).

F
I G U R E 4 (a) Boxplot showing the median (red line) Barlow ratio for the seven individuals.Boxes show interquartile range and whiskers show the lower and upper quartiles (Outliers are shown in Figure S2).(b) Bar plot showing the percentage of each portion of eye with a Barlow ratio that fell between 0.4 and 1.

F
I G U R E 5 (a) Cone shapes from all 7 individuals plotted by the length of the cone and the ratio of the diameter of the cone 10% from the bottom and the maximum width of the cone at the top.Scale bar shows 10 μm.Black line denotes the median cone length, while the red line shows the cone length of the typical cone.(b) Histogram showing the cone lengths across the eye of Polites peckius with example cones from the 3 different peaks in the distribution.(c) Histograms showing the aspect ratios and (d) cone ratios of each individual (excluding Heliconius as there were issues with individual cone segmentation in this dataset).